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O ' Abstract 
(N 

Ch , The distribution of masses for neutron stars is analyzed using the Bayesian statistical 

inference, evaluating the likelihood of proposed gaussian peaks by using fifty- four measured 
lO [ points obtained in a variety of systems. The results strongly suggest the existence of a bimodal 
' distribution of the masses, with the first peak around 1.37Mq, a much wider second peak at 



1.73M0. The results support earlier views related to the different evolutionary histories of the 
members for the first two peaks, which produces a natural separation (even if no attempt to 
"label" the systems has been made here), and argues against the single-mass scale viewpoint. 
Q^i The bimodal distribution can also accommodate the recent findings of ~ Mq masses quite 
Q ! naturally. Finally, we explore the existence of a subgroup around 1.25Mq, finding weak, if 
any, evidence for it. This recently claimed low-mass subgroup, possibly related to O—Mg—Ne 
j§ I core collapse events, has a monotonically decreasing likelihood and does not stand out clearly 
from the rest of the sample. 

I Keywords Stars: neutron 

00 ■ !• Introduction 

The measurements of masses and radii of neutron stars have the potential to constrain 
stellar evolution and dense matter physics alike. In fact, matching stellar evolution results 
and arguments to actual neutron star observations is crucial to test the whole theory. As is 
well-known, he standard Stellar Evolution theory suggests the mass at the main sequence to 
^ , be the most important parameter for the final outcome of massive stars. The lowest end of 
^1 ~ 8 — IIM0 is expected to produce very degenerate O — Mg — Ne cores which eventually 
collapse because of electron captures. Some possible systems in which the collapse can occur 
have been discussed in Siess (2007), Poelarends et al. (2008) , Nomoto &i Kondo (1991) 
Nomoto & Iben (1985) and Nomoto (1987), among other works. Podsiadlowski et al. (2005) 
elaborated on this problem and suggested that, due to the "characteristic mass" expected for 
the core, the equation of state and amount of ejected mass would produce low-mass neutron 
stars in the ballpark of ^ 1.25M0, also expected to receive small natal kicks in their birth 
events and therefore showing a low eccentricity. Further evidence in favor of this proposal 
has been presented by Schwab, Podsiadlowski & Rappaport (2010) after analyzing a sample 
of 14 well-measured neutron stars. In addition, a group of ~ 1.35M0 has been also identified 
and associated with the standard scenario of iron core collapse. This group features a much 
higher natal kick, and comprises some of the binary pulsar systems such as PSR 1913+16. 
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On the other hand, increasing evidence for massive neutron stars has been mounting, with 
several systems in the range ~ 1.6 — 1.8Mq and the very recent report of a 1.97 it O.O4M0 
(Demorest et al. 2010). Interestingly enough, most of these neutron stars typically have 
a white dwarf companion, and thus a distinct evolutionary history. However, systems in 
HMXRB and in a binary pulsar also exist. The mass could be the direct result of the 
existence of massive iron cores for M > 19Mq in the Main Sequence (Timmes, Woosley & 
Weaver 1996), a view advocated by van den Heuvel (2010) and others. 

The knowledge of this mass distribution is therefore fundamental to understand the mech- 
anisms involved in the final stages of stellar evolution. In this work, we consider the sample 
of know neutron star masses. We applied a Bayesian analysis for all set of masses overcome 
an a priori distribution. This work is motivated by a similar earlier approach by Finn (1994), 
who applied the Bayesian approach for to estimate the upper and lower limit for the neutron 
star masses, using one small set of data with only four well-measured neutron star binary 
pulsar systems. He observed the coincidence that all binary pulsar systems have constrained 
the masses close to 1.35M0. Schwab, Posialdowski & Rappaport (2010) recently argued for 
the existence of two distinct neutron stars populations, the first with high- mass ~ 1.35Mq 
and the second with low-mass ~ 1.25Mq, in a work that analyzed fourteen well-measured 
objects with uncertainties of < O.O25M0. They interpreted these two populations as to 
be result of distinct evolutionary formation scenarios: low-mass populations originates in 
electron-capture SNe and feature low kicks, while the high-mass population is the result of 
iron core collapse SNe. Both the situations, the authors proposed that masses are bounded 
the formation mechanisms restrict the range of neutron stars and the existence of the two 
channels for the production of neutron stars. Work by Thorsett and Chakrabarti (1999) a 
decade ago concluded that a very narrow range of masses around ~ 1.35Mq was found ana- 
lyzing a sample of 50 objects, ruling out accretion at the level of AM > O.IMq for the binary 
pulsars known at that time. Very recently, Zhang et al. (2010) concluded that a substantial 
accretion was present in recycled objects using an enlarged sample. Clearly, a reanalysis of 
this subject is in order. 

2. Neutron star sample 

Our adopted sample is the compilation by Lattimer and collaborators, publicly available 



at http : // st ellarcollapse . or g/nsmasses I . After finding the central values for each source (Fig. 
1), we included the 55 neutron stars with error bars varying between 0.009Mq and O.548M0 in 
our work. The full references to the original works, including the label letter employed in the 
compilation are listed in the References. The uncertainties for each mass are quite different 
because the methods of measurement were distinct at different times, and the methods had 
been improved in many cases. 

3. Statistical methodology 

In this work the analysis the sample of neutron star masses has been studied using 
Bayesian Statistics based on conditional probabilities, usually stated as p{Hi\D, I) = ; 
where p{Hi\I) is the probability priory p{D\Hi, I) is the likelihood function, p{Hi\D, I) is the 
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Figure 1: Neutron Stars masses sample. The corresponding uncertainties vary widely and 
have been included in the analysis. Data taken from stellarcollapse.org. 

posteriori probability and p{D\I) is the predictive probability. This approach has been shown 
to be particularly powerful for the treatment of scarce/inaccurate data, yielding nonetheless 
quite accurate estimate of the parameters in most cases. One of the main tools to evaluate the 
quality of the results is the Bayesian Information Criterion (BIG), which considers the ratio 
of the likelihood models, the available data and models, and explicitly penalizes the model 
having more parameters. A crucial (and often criticized) ingredient is the a priori expecta- 
tions from theoretical inference which is weighed by the BIG. An earlier study by Finn (1994) 
used radio observations for four neutron star binary pulsar systems and employed Bayesian 
Statistics, approximating each observed point by a gaussian function (on mass and standard 
deviation). He used one flat a priori distribution between an assumed upper limit of mass 
rriu and the lower limit mi. The values found were 1.01 < itii/Mq < 1.34 (lower limit) and 
1.43 < ttLu/Mq < 1.64 (upper limit). We aimed to improve this kind of analysis by working 
with the much larger sample and exploiting the potential of the Bayesian formalism. 

3.1 Likelihood 

Our tasks in this section is to construct the likelihood function. The likelihood distribu- 
tion is the key point of Bayesian analysis because it considers the data and the theoretical 
knowledge about the measurements together. Here, we assumed that the likelihood function 
is simply the product of independent probabilities for what was measured and what was 
expected to be measured 



r 

L{9\D,M) = 



l[pim,\D,M) 



l[pim,\D,M). 



Where 9 represents the space of parameters, D is the data set, M the models, p{mi\D, M) 
is the probability of the data to be measured and p{mj\D, M) is the expected probability for 
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the measured data. The hkehhood weights the samphng probabihty given the data (D) and 
assuming the model (M). The hkehhood in our case is 

L{e\D,M) = - exp J np{M,Mi,M2,ap,ai,a2,G)dM 

YlapX g{Mi,ai,Mi,Ci) + (1 - ap - ao)x 

xg{M2,a2,Mi,i,) + x 5(1.25, 0.07, Mi , ^i) 

Where Up is a function that involves the peak masses Mi, M2, ap is the relative amplitude 
of the first peak, ac0 is the amplitude centered on the 1.25M0, with assumed standard 
deviation ~ 0.07, (Ti, o"2 are the standard deviations of the theoretical peaks and g \s a. 
gaussian function (actually the product of two gaussian distribution integrated). 



9=1 exp 

/ ai 



{u — XiY 



X exp 



(n - X2)^ 



2ql 



du 



Where ai and 02 are 



and 



oi = max[xi — H X qi,X2 — H X (72] 



02 = max[xi + H X qi, X2 + H X q2]; 
H is the scale parameter, qi , q2 are the standard deviations of the xi and X2- Finally, 



np{M,Mi,M2,ap,ai,a2,g) = ap x exp 
(M- 1.25)2^ 



{M -Mif 



2si 



+ 



+G X exp 



2(0.072) 



+ [1 - (op + ao)] X exp 



(M - M2f 
2al 



3.2 The a priori distributions 

The a priori distribution is the assumed knowledge about the phenomenon that is treated. 
In many previous works, notably those of Finn (1994) and Schwab, Podsiadwolski and Rap- 
paport (2010), gaussian distributions were employed, and are a "natural" choice employed in 
our work (see Finn 1994 for a discussion of this gaussian form within the Bayesian approach). 
We then assumed, to make contact with the work of Schwab, Podsiadwolski and Rappaport 
(2010), a distribution peaked on two masses, around Mi = 1.35M0, and M2 = 1.55AfQ. 
Those authors restricted their analysis to a set of 14 well-measured NS, and did not include 
the large error bar points of the X-ray binaries and several WD-NS systems, consistently 
with their frequentist approach. With the aim of exploring the full distribution, we have 
chosen a second, higher value of M2 to match the plain mean value of the NS in WD-NS 

^Where < Op + ao < 0. 



4 



and X-ray binary systems of the sample first but, as we shall see below, the precise value is 
not too important at this point. We performed a first run within this two-value hypothesis, 
and compared it with the occurrence of a single mass peak for all objects, including NS-NS 
binaries, WD-NS binaries, X-ray binaries and MS-NS systems. The motivation for this first 
bold comparison was to see whether a single mass scale was still possible with the present 
data, as many works insisted on till a few years ago. 

After this "first run", and still working within the two-peak hypothesis, and adopting the 
same values as Schwab, Podsiadwolski and Rappaport (2010) for Mi and o"i, we looked for 
evidence of a subgroup attributed by them to O — Mg — Ne cores producing low- mass neutron 
stars. In addition to the main peaks, we examined the distribution to identify the presence 
of this subgroup by defining oq as the peak amplitude of a 1.25Mq, that is, relative to the 
two main identified peaks. Again, since the claimed masses lie at the low end, the calculation 
is not sensitive to the precise value of M2. We defined and calculated the relevant likelihood 
of a third peak, physically related to the formation of light NS out of O — Mg — Ne cores of 
the lighter progenitors, of this second run. 

In a third run we left the values of the masses Mi and M2 to be determined directly by 
the raw available data, together with their standard deviations values ai and a2- That is, 
while still imposing the distribution to be composed of two gaussian forms, we sought for 
the optimal values without restrictions as driven by the data sample. The purpose here is to 
let the Bayesian tools to indicate which are the possible values given the full error bars and 
within a definite gaussian hypothesis. 

4. Results and Conclusions 

We first address the basic results provided by the first run: the likelihood of still a single 
mass-scale as an explanation of all the data points is much lower than (at least) two peaks in 
them, in spite of the introduction of extra parameters penalized by the Bayesian approach. 
This results may not be meaningful for some, since it has been known for years that the more 
massive systems should have accreted 0.1 — O.2M0 or more, and therefore detach from the 
original mass value. However, the first run is important to overcome the idea that a single 
value of the NS mass will be enough: even if the two peaks are quite close, this distribution 
is preferred to a single wide peak. This is somewhat expected, since the extremes of the 
determined masses, around IMq, (van der Meer et al. 2007, see the recent work by Rawls 
et al. 2011 released after the completion of this work) at the lowest and 1.97 it 0.04Mq 
(Demorest et al. 2010), with a few more massive candidates (Clark et al. 2002, Preire et al. 
2008b), are now separated by at least IMq. 

The results of the second run are summarized in Fig. 2. The calculated likelihood of the 
third peak is a monotonically falling function of the amplitude, being maximal around zero 
amplitude (relative to the prominent main peaks). Therefore, we conclude that there is little 
evidence for the presence of a low-mass subgroup, and that the four objects in this range may 
be in fact members of the 1.35M0 peak . However, the confirmation of this result would have 
important consequences, since the progenitor stars are quite abundant. The lack of strong 
evidence of a peak ~ 1.25M0 could be due to still poor sample/bias, or it could alternatively 
mean that most of the 8 — IIM0 become AGB, an issue that deserves serious consideration. 
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Figure 2: The likelihood of a third, low-mass peak. The existence of this third peak has a 
decreasing likelihood except for a nearly constant value near oq ~ 

The results of the third run rendered two masses, as shown in Fig. 3. The first mass is not 
substantially different from the result of Schwab, Podsiadwolski and Rappaport (2010), a fact 
we interpret as the robustness of the distribution at this scale; while the second one is now 
around 1.73Mq. This is not surprising, since the known WD-NS and X-ray binary systems 
have large uncertainties but also high values of the central NS masses. The full advantage of 
the Bayesian techniques suggest here that masses around 2Mq are not unexpected, even less 
when considering the obtained value of (T2 = 0.25 (for comparison, the value of ai = 0.042 
is very narrow, as expected). The recently announced value of the object in the system, 
1.97 ± 0.04Mq according to Demorest et al. 2010 (not included in the sample), is just an 
example of this higher figure for the systems of this kind, and suggests a mean accreted 
mass of several tenths of Mq, although the precise value is model-dependent and should be 
estimated for each individual case. This subgroup may also include members coming from 
> 2OM0 masses in the Main Sequence (Timmes, Woosley and Weaver 1996, van den Heuvel 
2004), born with higher masses without suffering any substantial accretion.. 

In summary, we have presented an analysis of an available NS mass sample which indicates 
a) a bimodal distribution with a narrow peak at Mi = 1.37Mq (fully compatible with the 
findings of Schwab, Podsiadwolski and Rappaport (2010), and a higher peak at Mi = I.TSMq, 
the latter with a wide shape capable of accommodating > 2Mq masses, b) little evidence for 
the expected low-mass NS descendants from the 8 — IIMq range in the Main Sequence, and 
c) the unsustainability of the "one-mass-fits-all" picture, since the adopted Bayesian scheme 
properly weights the large uncertainties and extra parameters of the bimodal hypothesis 
applied to an imperfect dataset. Given the potential implications for the stellar evolution, the 
need of further pursuing these kind of analysis and a continuous improvement by incorporating 
new data can not be overstated. 

After the completion of this work we noticed the public release of a Bayesian analysis by 
Kiziltan, Kottas and Thorsett (2010), tackling the same problem and with a very detailed 
treatment of the calibration. Also their sample is restricted to avoid the inclusion of WD- 
NS and X-ray binary systems with large error bars. Within this analysis, the authors find 
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Figure 3: The values of masses at which the distribution peaks. The numbers on the contours 
indicate the decreasing probabilities of finding Mi, M2 there. 



two peaks, one at M = I.SSMq, and a second one located at M = 1.5Mq. Given the 
different sampling and their reliance on simulated distributions (not performed by us), it is 
fair to say that there is enough room for a convergence of results. It is important to note 
that these authors warn against the inclusion of more uncertain points because of a possible 
contamination of the sample. While we agree with this judgement, we believe that unless the 
higher masses happen to be completely and systematically overestimated, the emergence of 
a bump at a mass even higher than I.SMq is strongly expected. Our value for that "second 
peak" should reflect to a large extent the difference of including these loose observational 
data as processed by the Bayesian formalism, acknowledged to handle this kind of situations 
better than well-known frequentist approaches. In addition, the quest for a low-mass bump 
remains in that work, since the peak at the lower M = 1.35Mq hosts the low-mass systems 
attributed by Schwab, Podsiadwolski and Rappaport (2010) to the lower end of progenitor 
masses. 

Finally, we acknowledge the works of Steiner, Lattimer and Brown (2010) and Zhang et 
al. (2010), the first aimed to reveal the nuclear equation of state and the second focusing on 
the evolutionary features of the systems. Steiner, Lattimer and Brown found independently, 
from an analysis of a subset of X-ray bursters, that the maximum mass had to be quite high, 
as demanded by the 2-Mq determination by Demorest et al. (2010) released shortly after 
their paper. We did not attempt any specific inference about the nuclear equation of state 
here, although our results also demand a theoretical description capable to accommodate the 
objects in the second peak. On the other hand, Zhang et al. used almost the same sample, 
but their calculations attempted to link the period to the mass, rather than discriminating 
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between single-peak and multimodal distributions. Even their statements about the NS-NS 
systems alone does not allow a firm conclusion about the existence of the low-mass peak at 
1.25Mq, since they focus just on the mean values and dispersions. In all cases the relatively 
high values of the latter dispersions seem to be a consequence of this methodology. 
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